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The time dependence of the survival probability, S(t), is determined for diffusing particles in 
two dimensions which are also driven by a random unidirectional zero- mean velocity field, v x (y). 
For a semi-infinite system with unbounded y and x > 0, and with particle absorption at x = 0, a 
qualitative argument is presented which indicates that S(t) ~ t _1//4 . This prediction is supported 
by numerical simulations. A heuristic argument is also given which suggests that the longitudinal 



probability distribution of the surviving particles has the scaling form P(x,t) ~ t 



d/3 



g{u). Here 



the scaling variable u oc x/t 3 ^ 4 , so that the overall time dependence of P(x,t) is proportional to 
t~ B ^ 4 , and the scaling function g(u) has the limiting dependences g(u) oc const, as u — * and 
g(u) ~ exp(— M 4,/3 ) as u — > oo. This argument also suggests an effective continuum equation of 
motion for the infinite system which reproduces the correct asymptotic longitudinal probability 
distribution. 



PACS Numbers: 05.40. +j, 05.60.+W, 02.50.Ey 



I. INTRODUCTION 



Consider a diffusing particle in the semi-infinite pla- 
nar domain (x > 0,y) which is absorbed when x = 
is reached. In addition to the diffusion, the particle is 
driven by a unidirectional random velocity field in which 
v x (y) is a random, zero-mean function of y only (Fig. 1). 
This type of stochastic motion was introduced by Math- 
eron and de Marsily (MdM) |Q to describe the hydrody- 
namic dispersion of dynamically-neutral tracer in a sed- 
imentary layered rock formation. Although the longi- 
tudinal bias averaged over an infinite number of trans- 
verse layers is zero, the typical bias over a finite number 
of layers is a fluctuating quantity which is a decreasing 
function of the number of layers that a random walk vis- 
its. This non-vanishing residual bias underlies the faster 
than diffusive transport of the model. In an infinite sys- 
tem it has been established that the typical horizontal 
displacement x typ oc £ 3 / 4 ||-|^] and that the configura- 
tion averaged distribution of longitudinal displacements 
has the form P{x,t) ~ r 3 / 4 exp[-(x/i 3 / 4 ) 4 / 3 ]. There 
are also strong fluctuations in the probability distribu- 
tion between different samples of the velocity field, as 
well as a very slow convergence to the asymptotic limit. 

While much is now understood about transport in the 
MdM model, we wish to investigate its first passage prop- 
erties. Specifically, we consider the semi-infinite two- 
dimensional system x > with an absorbing boundary 
at x = 0, and study the time dependence of the parti- 
cle survival probability, S(t). The survival probability 
in a finite system with absorbing boundaries at x = ±L 
and the same unidirectional random velocity field has 
been studied previously Q ; however this system exhibits 
fundamentally different behavior than the semi-infinite 
system that is treated here. 

In the absence of a velocity field, it is well known that 
in the semi- infinite system S (t) asymptotically decays in 
time as i~ 4 / 2 |5). Because the velocity field in the MdM 



model has no net longitudinal bias, it is not immedi- 
ately obvious how the behavior of S(t) will be affected. 
Naively, one might expect that the dominant contribu- 
tion to S(t) will arise from those velocity configurations 
whose average bias is directed away from the boundary. 
This is indeed the case, and from this starting point, 
we present a simple argument which suggests that the 
survival probability, averaged over all realization of the 
velocity field, is proportional to t -1 / 4 . This prediction is 
in excellent agreement with our numerical results. 



x=0 



FIG. 1. The random velocity field in a realization of the 
MdM model on a finite width square lattice strip. On the 
horizontal bonds, the direction of the velocity field is indi- 
cated by the arrows. In a single time step, a particle (heavy 
dot) can move equiprobably only to one of the two target sites 
indicated. Particles are absorbed at x — 0. 

Interestingly, a i -1 / 4 decay of the survival probabil- 
ity has been found previously for diffusing particles in 
a semi-infinite two dimensional system with a unidirec- 
tional zero-bias, but deterministic velocity field of the 
form v x (y) = —v x (—y) Although the mechanism 

that leads to S(t) ~ i" 1 ' 4 in this class of velocity fields is 
different than that for the MdM model, the two systems 
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share the feature that their velocity fields have no average 
bias in x. It would be interesting to determine whether 
a i -1 / 4 decay of the survival probability is characteristic 
of all semi-infinite systems with a zero-mean longitudinal 
velocity field. 

In the next section, we give a heuristic argument which 
suggests that S(t) oc i -1 / 4 . We then present corroborat- 
ing numerical simulations in Sec. III. We also find that 
the spatial probability distribution of the surviving parti- 
cles provides insight into the formulation of a continuum 
equation of motion for the longitudinal probability dis- 
tribution in an unbounded geometry. Thus in Sec. IV, we 
infer this equation of motion and, from its solution, de- 
termine the correct asymptotic longitudinal probability 
distribution in the unbounded geometry. We conclude 
with a brief discussion in Sec. V. 



II. CONFIGURATION-AVERAGED SURVIVAL 
PROBABILITY 

We first present our argument for the the time de- 
pendence of S(t). The basic idea is that the dominant 
contribution to this average arises from the subset of all 
velocity configurations whose net bias is away from the 
boundary, i. e., in the +a;-direction. Conversely, config- 
urations with a bias along —x will give individual con- 
tributions to S(t) which decay exponentially in time and 
thus should be asymptotically negligible. 

To determine which of the positively-biased velocity 
configurations give the dominant contribution to S(t), 
consider a discrete realization of the MdM model on the 
square lattice in which the velocity is either +vq or — vq 
with equal probability for a given value of y (Fig. 1). Pe- 
riodic boundary conditions in the transverse direction are 
employed, so that the system is a semi-infinite cylinder 
consisting of w rows. For concreteness, the initial con- 
dition is p(x,y,t = 0) = ^3x,e, «■ e., a ring of particles 
is initially placed at x = £, where £ is the lattice spac- 
ing. In this system, the probability that there are n+ 
positively biased rows and n_ negatively biased rows is 
V(m) cx -y=e~ m l w , with m = n + — n_. In a time t, the 

number of layers visited by a random walk is w oc \f~Di/i, 
where D is the transverse (microscopic) diffusivity. By 
transforming from m to the velocity v — rnvg/w, the 
distribution of velocities within w layers is 
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Because this distribution is strongly cut off when the ar- 
gument of the exponential is greater than one, we expect 
that the dominant contribution to S(t) will arise from 
those velocity configurations whose net bias is within the 
range < v < v (£ 2 / Dt) 1 / 4 . 

For a positively biased velocity configuration, we now 
estimate the residual survival probability in the long time 



limit under the assumption that w is finite. In this case, 
the particle will uniformly sample the transverse extent 
of the system and it is sensible to characterize the bias 
by its mean value v. If v is small, or more properly 
the Peclet number vxq/D is small, then for t < D/v 2 
the bias is irrelevant and consequently S(t) ~ x$j\J Dt, 
where xq is the initial position of the particle |q]. How- 
ever, for t > D/v 2 , convection dominates and returns 
to the origin become extremely unlikely. Hence S(t) 
should "stick" at the value attained when t = D/v 2 . 
This implies that the asymptotic behavior of the survival 
probability is simply, S(t = oo) oc vxq/D. This same 
result can also be obtained with additional effort from 
a more rigorous approach in which one solves the one- 
dimensional convection-diffusion equation on the domain 
x > with the initial condition P(x,t — 0) = 5{x — xq) 
and then computes S(t) by integrating this probability 
density over all x > 0. 

For the MdM model, we now average over all relevant 
velocity configurations to determine S(t). That is 
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III. NUMERICAL SIMULATIONS 

We now present numerical evidence to support the pre- 
diction that S(t) oc i -1 / 4 . For simulating the particle mo- 
tion, we propagate the probability distribution for each 
velocity configuration exactly 0. The microscopic mo- 
tion is defined by the rule that in a single time step, a 
particle hops equiprobably by ±1 in the y-direction and 
by a distance sign(v(y))£ in the x-direction on the square 
lattice, where v(y) is the velocity field at the y-coordinatc 
of the particle before the hopping event (Fig. 1). Accord- 
ingly, the probability that a particle is at [x, y) at time t 
evolves according to 

p(x, y,t+l) =ip(x-£x sign(«(y -£)),y-£, t) 

+ \p{x-£x sign(«(y + £)),y + £, t). (3) 

For each velocity configuration, probability propagation 
yields the exact distribution in the presence of the ab- 
sorbing boundary up to the maximum time specified. 

We have performed the average over velocity configu- 
rations in two complementary ways - either an average 
over a representative sample of velocity configurations, 
or an exhaustive average over all velocity configurations 
for relatively small systems. The former is straightfor- 
ward to implement, but the effect of rare configurations 
on the results is unknown. The latter, on the other hand, 
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gives the exact result for an infinite system, albeit only 
for short times. This exactness allows one to test for 
systematic trends in the data, an analysis which is not 
feasible when averaging over a representative set of ve- 
locity configurations. 

A typical result for S(t) up to t = 4095, based on 
an average over 50 velocity configurations on a cylinder 
of width 400, is shown in Fig. 2. The average bias of 
these 50 configurations turns out to be —0.004. Beyond 
approximately 50 time steps, the data for the survival 
probability is quite linear and a least-squares power law 
fit to the data in this time range yields the exponent 
of —0.2491. Further, the slope between successive data 
points, or local exponent estimate 



\n{S{t)/S{t-l)) 
w[) ~ In(t/(t-l)) • 



(4) 



deviates from —0.25 by less than 0.002 for t > 50 (in- 
set). To check that finite width effects do not substan- 
tially affect the results, we also considered shorter times, 
t < 1600 and a slightly narrower system (w = 300), 
and averaged over 250 realizations. The average bias of 
these configurations turns out to be +0.002. This case 
yielded a best fit exponent of —0.2503. These two data 
sets strongly suggest that S(t) oc t -1 / 4 in the long time 
limit. However, because the average is performed only 
over an infinitesimal fraction of all velocity configura- 
tions, it is possible that extreme configurations could al- 
ter the results. For this reason, we now investigate the 
exact behavior of S(t) for short times by averaging over 
all velocity configurations. 
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In this complete enumeration, we consider odd values 
of w for convenience. Because of the periodic transverse 
boundary conditions, many of the 2 W configurations are 
identical up to cyclic permutation and reflection symme- 
try. To carry out the enumeration, we first encode each 
velocity configuration as a binary sequence. Using bit 
manipulation techniques, we identify the "irreducible" 
representation of this binary sequence, defined as the 
smallest equivalent integer number obtained by perform- 
ing all possible cyclic permutations of the initial binary 
sequence. This same procedure is then repeated on the 
reversed initial binary sequence. Thus to each binary 
sequence there is a unique irreducible binary sequence. 
By this mapping, we only need consider the irreducible 
configurations and weight each by their degeneracy in 
performing the average over velocity configurations. For 
example for w — 23, 25, 27, and 29, the number of irre- 
ducible configurations are 92,205, 337,594, 1,246863, and 
4,636,390, compared, e. g., to 2 29 = 536, 8 70, 9 1 2. For 
each irreducible configuration, we then perform the ex- 
act probability distribution propagation. This complete 
enumeration provides the exact value of S(t) for an in- 
finite system up to w — 1 time steps, while finite width 
crossover effects gradually begin to play a role for later 
times. 
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FIG. 3. The local slope of S(t) versus ! on a double loga- 
rithmic scale, based on averaging the exact probability distri- 
bution over all of velocity configurations for systems of width 
w = 15, 17, . . ., 29. The data for different w coincide for 
larger 1/t, but then separate at progressively smaller 1/t as 
w increases 



FIG. 2. Plot of S(t) versus t on a double logarithmic scale, 
based on averaging the exact probability distribution over a 
finite number of velocity configurations for the system dis- 
cussed in the text. The inset shows the local slopes between 
neighboring data points. 



We therefore typically carried out the probability prop- 
agation for up to t w 2w time steps and exploited the 
crossover in S(t) to interpret our results. A plot of In S(t) 
versus lnt should initially show power law behavior, in- 
dicative of the infinite system behavior, and then cross 



3 



over to a non-zero constant because of finite width ef- 
fects. Thus a plot of the local exponent a w (t) (here de- 
fined as the slope between every other data point) versus 
1 jt should initially provide an estimate of the exponent of 
S(t), while the crossover effect determines the time range 
over which the exact data is relevant for the infinite sys- 
tem. In Fig. 3, this local slope, is plotted versus 1/t for 
system widths w between 15 and 29. Initially, a w (t) is 
decreasing nearly linearly in 1/t, but subsequently there 
is the expected crossover to the asymptotic value of zero. 
In the regime where the data is relatively linear, we com- 
pute the intercepts of successive data points at 1/t = 
as an estimate of the asymptotic value of the exponent 
(Fig. 4). As w increases, this data exhibits: (i) non- 
monotonic trends in the data (e. g., the location of the 
minimum) which disappear only for w > 25, (ii) more 
stable extrapolated values as w increases, and (iii) the 
minimum value of the extrapolated exponent - which we 
adopt as the best estimate of the exponent for a given 
value of w - is increasing systematically in w and appears 
to be converging to —0.25. While there is slow conver- 
gence to asymptotic behavior which gives rise to consid- 
erable subjectivity in analysis, we believe that the trends 
in the data support the hypothesis that S(t) ~ i -1 / 4 . 
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FIG. 4. Linear extrapolation of the local slope from Fig. 3 
for w — 15, 17, . . ., 29. The minimum value is progressively 
increasing with w. 



In addition to investigating S(t), we also examined the 
probability distribution of the surviving particles. This 
quantity provides an alternative understanding for the 
first-passage process, as well as useful fundamental in- 
sights about the continuum description of MdM model. 
Specifically, we study P(x,t) = J p(x,y,t) dy, the con- 



figuration averaged longitudinal probability distribution 
of particles which have not yet been absorbed by time 
t. We expect that this probability distribution can be 
written in a scaling form 

P(x,t) = Af(x/(x)), (5) 

where (x) — (x(t)} is the average longitudinal displace- 
ment of the survivors at time t. Monte Carlo simula- 
tions clearly indicate that (x) cx as in the case when 
there is no absorbing boundary present Because 
J P(x, t) dx = S(tY. we can determine the coefficient A 
by integrating Eq. (|5J) over x and thereby write 

P(x,t) = ^rf(x/{x)), (6) 

where Fq = L dufiu). Because of the absorbing 
boundary, P(x = 0, t) must equal zero, leading to the ex- 
pectation that f(u) will vanish as a power law as u — > 0. 
Consequently, we write 

p(x ' t)= (KM)^ W(i>) ' (7) 

with g{u) — > const, as u — > 0, and g(u) vanishing faster 
than any power law for u — > oo. 

A plot of the scaling function f(u) versus u is shown 
in Fig. 4 for t = 255, 1023, and 4095. There is a small 
but systematic variation in the data for different times, 
with the small-w behavior steepening and the large-u tail 
growing for larger time. Nevertheless, reasonable data 
collapse is obtained in which f(u) qualitatively exhibits 
the expected power-law and rapid cutoff asymptotic be- 
haviors for small and large u, respectively. We attribute 
the small deviation from scaling on slow convergence to 
the asymptotic limit (see below). Such a phenomenon 
was observed previously in the probability distribution 
for an infinite system [J^H , and similar slow convergence 
effects can be anticipated here as well. 
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IV. THE PROBABILITY DISTRIBUTION 0.5 



4 



FIG. 5. The scaling function f(u) versus u. 

The exponent /i in Eq. (Q) can be obtained by demand- 
ing consistency between the time dependence of S(t) and 
that of the first-passage probability. If S(t) oc i -1 / 4 , 
then from the general relation ^ between S(t) and the 
first passage probability to the boundary, T(t), we have 
^ fi^ ■ On the other hand for a normalized initial 
condition, the first passage probability coincides with the 
flux to x = 0. Now as x — > 0, Eq. (Q) gives 



n , \ S(t) 
P(x,t) oc ar oc 



t* ^ £-1/4-3(1+/*) /^M. 



(8) 



Since the flux is obtained by performing an appropriate 
spatial derivative of this limiting probability distribution, 
an operation which does not affect the temporal behav- 
ior, we conclude that fx = 1/3 to recover the correct t~ 5 / 4 
time dependence for the flux. 

However, the data in Fig. 5 does not exhibit this be- 
havior because of finite time effects. For small u, the 
numerical value of g(u) at x = 1 is non-zero but de- 
creasing with time. Correspondingly, the value of u at 
this first data point is non-zero but also decreasing with 
time. This anomaly in the small-u data renders a simple 
power law fit inadequate. However, such a naive fit to 
the data in the range u < 1/2 gives the estimates 0.483, 
0.479, 0.476, 0.470, and 0.464, respectively, for the expo- 
nent n in P(x, t) for the 5 aforementioned time values. 
Other analyses, such as computing the first derivative of 
/(it) (which should diverge as u~ 2 / 3 ) and determining 
the exponent of the small-it dependence, lead to a simi- 
lar quantitative conclusion. While the anticipated value 
/i = 1/3 is not obtained, we believe that better agreement 
with theory would emerge if it were practical to extend 
our Monte Carlo simulations to much longer times. 

The existence of this power-law prefactor further sug- 
gests that in a continuum description, the configuration 
averaged flux can be obtained via 
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r 2/3 



dP(x, t) 



dx 



(9) 



Such a scale-dependent diffusivity can be heuristically 
justified by appealing to a Taylor diffusion description of 
the motion in the MdM model B. In a time t, a par- 
ticle typically explores w « \f~Dt transverse layers. The 
typical bias within this number of layers is then propor- 
tional to w~ x l 2 or (Z?t) -1 / 4 . Thus in a time scale t, 
the typical longitudinal distance travelled by a particle 
is d ~ vt oc i 3 / 4 . Because these segments of length d 
are randomly in the +x or — x direction and the time 
interval between segments is of order t, we infer an ef- 
fective longitudinal diffusion (or dispersion) coefficient 
D|| - £ 2 /t oc t 1 / 2 oc x 2 / 3 , as written in Eq. (§). 

Let us now pursue the consequences of this scale de- 
pendent diffusion coefficient for the longitudinal motion. 
If the longitudinal flux j(x,t) is indeed proportional to 



—x 2 ^ 3 dP(x t) • • • 

Q x y ' ' , then substituting this into the continuity 

equation 



dP(x,t) , dj(x,t) 



0. 



dt dx 
leads to the effective equation of motion 

/3 dP(x,t) 



dP( X ,t) = d_ x 2/3 



dt 
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(10) 



(11) 



We can easily solve this equation by applying scaling. 
Assuming that P(x,t) oc t~ 3 / 4 h(x/t 3 / 4 ), we rewrite the 
partial derivatives in x and t in terms of a derivative with 
respect to u = x/t 3 / 4 , to recast the equation of motion 
as 



'-(uh(u))' = (u 2/3 h'{u))'. 



(12) 



Here the prime denotes differentiation with respect to it. 
One integration immediately yields 



-uh{u) = (u 2/3 h'(u)). 



(13) 



The constant of integration equals zero because h{u) — > 
faster than any power law as it — -> oo. A second integra- 
tion then gives h(u) oc exp(— it 4 / 3 ), from which we con- 
clude that the longitudinal probability distribution has 
the form 



P(x,t) oc i- 3 / 4 exp(-(a;/< 3/4 ) 4/3 ). 



(14) 



This functional form coincides with that obtained previ- 
ously by a different method |2|,|3| in which the dominant 
contribution to the large-u tail of P{x 1 1) arises from ex- 
treme "stretched" trajectories in unlikely velocity con- 
figurations. Thus the observation of the large-u tail for 
P(x, t) from the numerical data in Fig. 5 can again be an- 
ticipated to be problematical; much more extensive sim- 
ulation would be needed. 



V. DISCUSSION 

We have investigated the time dependence of the con- 
figuration averaged survival probability, S(t), in a semi- 
infinite two-dimensional system for diffusing particles 
which are also driven by a unidirectional random zero- 
mean velocity field, v x (y). A qualitative argument sug- 
gests that S(t) oc i -1 / 4 , a prediction which is in excel- 
lent agreement with numerical results. We also examined 
the longitudinal spatial probability density of the surviv- 
ing particles, P(x,t). Interestingly, although the numer- 
ical evidence supporting the prediction that S(t) oc i~ 4 / 4 
is strong, the numerical data for P(cc, t) indicates slow 
convergence to the scaling limit and some inconsistency 
with the behavior of S(t) itself. Similar anomalies in 
the probability distribution occur in the unbounded ge- 
ometry (|-||] , due to the contribution of extreme velocity 
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configurations in the average. It is surprising that the be- 
havior of S(t) is apparently relatively insensitive to the 
contribution of such extreme configurations. 

An interesting byproduct of the spatial distribution of 
the surviving particles is that the form -^|r^r^ is sug- 
gested as the appropriate expression for the particle flux. 
This leads to a scale dependent diffusion coefficient which 
is proportion to x 2 ^ , as well as the continuum equation 
of motion, Eq. ([ll]), for the longitudinal spatial probabil- 
ity distribution in an unbounded geometry. The solution 
to this equation of motion is simple to obtain and repro- 
duces the known asymptotic form of P(x, t) in the un- 
bounded geometry [|]-[| ■ As an application of this equa- 
tion of motion, we find new predictions for steady-state 
transport properties. For example, for a steady input of 
particles at x = L and particle absorption at x = 0, the 
steady solution to Eq. ([ll]) gives a configuration-averaged 
density profile which varies as X 1 / 3 . It will be worthwhile 
to test this prediction and also the general prescription 
for obtaining an effective equation of motion from the be- 
havior of the particle flux near an absorbing boundary. 

Finally, the behavior of S(t) for a semi-infinite system 
with a longitudinal MdM velocity field can be easily gen- 
eralized to arbitrary spatial dimension d. From classical 
results [pi , the number of distinct longitudinal rows vis- 
ited by a random walk in time t varies as t( d_1 '/ 2 f or 
dimension 2 < d < 3 (i. e., transverse spatial dimension 
between 1 and 2), as i/lni for d = 3, and as t for d > 3. 
Following closely the approach in Sec. II, this then leads 
to 

-1/4 d = 9, 



S(t) 



r 

t -{d-i)/i 2<d<3; 
{Xnt/t) 1 ' 2 d=3; 

d>3. 



(15) 



Thus above three dimensions, the survival probability ex- 
ponent value is not affected by the presence of a random 
velocity field. 
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